#####看727个位点的基因能否富集在gene-set数据库中，用基因名富集

##结果：727个位点，共822个基因，可以在FMRP中得到富集

library(meta)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

library(openxlsx)
fn=list.files(pattern = "xlsx",path="E:/0 公共数据库差异情况/db_for_5hmc/geneset",full.names=T)

tmp=read.xlsx(fn[1],sheet=1)
FMRP=data.frame(symbol=tmp$Gene.Symbol,group=rep("FMRP",842))

tmp=read.xlsx(fn[2],sheet=1)
GABA=data.frame(symbol=tmp$gene.symbol)
GABA$group="GABA"

tmp=read.xlsx(fn[3],sheet=1)
NMDAR=tmp[!is.na(tmp$NMDAR),]
NMDAR=data.frame(symbol=NMDAR$gene,group=rep("NMDAR",61))
ARC=tmp[!is.na(tmp$ARC),]
ARC=data.frame(symbol=ARC$gene,group=rep("ARC",28))


tmp=read.xlsx(fn[4],sheet=1)
tmp=tmp[!tmp$Approved.Gene.Name=="NOT_FOUND",]
PSD=data.frame(symbol=tmp$Approved.Gene.Name,group=rep("PSD",1447))

tmp=rbind(FMRP,GABA)
tmp=rbind(tmp,NMDAR)
tmp=rbind(tmp,ARC)
tmp=rbind(tmp,PSD)

 group=unique(tmp$group)
 

case=file2
caseid=unique(unlist(strsplit(as.character(case$Gene.refGene),";")))
caseid=caseid[!caseid=="NONE"]


#case1=merge(case,sig_deg,by.x = "Gene.refGene",by.y="symbol")


con.file=list.files(pattern="con_genotype.random",path="./20201207/",full.names=T)

col_names=c("geneset","con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
for(j in group){
tmp1=tmp[tmp$group==j,]
tmp1id=tmp1$symbol
for(i in 1:3){
a=length(intersect(tmp1id,caseid))
b=length(caseid)-a
con=read.table(con.file[i],head=T,sep=",")
conid=unique(unlist(strsplit(as.character(con$Gene.refGene),";")))
conid=conid[!conid=="NONE"]
c=length(intersect(tmp1id,conid))
d=length(unique(conid))-c
rt_tmp=data.frame(matrix(NA,1,ncol=8))
names(rt_tmp)=col_names
rt_tmp[,1]=j
rt_tmp[,2]=gsub("./20201207/","",gsub(".hg19_multianno.csv","",con.file[i]))
rt_tmp[,3]=a
rt_tmp[,4]=b
rt_tmp[,5]=c
rt_tmp[,6]=d
rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
rt_tmp[,8]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,rt_tmp)
}
}
 write.csv(result[-1,],"./20201231.GeneSet/GeneSet_enrichment.807psyASH.csv",quote=F,row.names = F)

 ###马德，meta分析不能用于同一个研究（case不能相同，以下作废，保存以上3种结果。）
result=result[-1,]
col2_names=c("geneset","con.file","case_overlap","meta.OR","meta.p.value","lower","upper")
rt.meta=data.frame(matrix(NA,1,ncol=7))
names(rt.meta)=col2_names
rt.meta=rt.meta[-1,]

for(i in group){
ftmp=result[result$geneset==i,]
ftmp$case_all=ftmp$case_overlap+ftmp$case_not
ftmp$con_all=ftmp$con_overlap+ftmp$con_not
rt.meta.tmp=data.frame(matrix(NA,1,ncol=7))
names(rt.meta.tmp)=col2_names
rt.meta.tmp[1,1]=i
rt.meta.tmp[1,2]="ramdon.150K"
rt.meta.tmp[1,3]=ftmp[1,"case_overlap"]
metaor=metabin(case_overlap,case_all,con_overlap,con_all,data=ftmp,sm="OR",studlab = rep(i,3))
rt.meta.tmp[1,4]=exp(metaor$TE.fixed)
rt.meta.tmp[1,5]=metaor$pval.fixed
rt.meta.tmp[1,6]=exp(metaor$lower.fixed)
rt.meta.tmp[1,7]=exp(metaor$upper.fixed)
rt.meta=rbind(rt.meta,rt.meta.tmp)
}